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Abstract — Recent work by Zymnis et a/, proposes an efficient 
primal-dual interior-point metliod, using a truncated Newton 
metliod, for solving tlie network utility maximization (NUM) 
problem. This method has shown superior performance relative 
to the traditional dual-decomposition approach. Other recent 
work by Bickson et al. shows how to compute efficiently and 
distributively the Newton step, which is the main computational 
bottleneck of the Newton method, utilizing the Gaussian belief 
propagation algorithm. 

In the current work, we combine both approaches to create 
an efficient distributed algorithm for solving the NUM problem. 
Unlike the work of Zymnis, which uses a centralized approach, 
our new algorithm is easily distributed. Using an empirical 
evaluation we show that our new method outperforms previous 
approaches, including the truncated Newton method and dual- 
decomposition methods. As an additional contribution, this is the 
first work that evaluates the performance of the Gaussian belief 
propagation algorithm vs. the preconditioned conjugate gradient 
method, for a large scale problem. 

I. Introduction 

We consider a network that supports a set of flows, each of 
which has a nonnegative flow rate, and an associated utility 
function. Each flow passes over a route, which is a subset 
of the edges of the network. Each edge has a given capacity, 
which is the maximum total traffic (the sum of the flow rates 
through it) it can support. The network utility maximization 
(NUM) problem is to choose the flow rates to maximize the 
total utility, while respecting the edge capacity constraints 
[1], [2]. We consider the case where all utility functions 
are concave, in which case the NUM problem is a convex 
optimization problem. 

A standard technique for solving NUM problems is based 
on dual decomposition [3], [4]. This approach yields fully 
decentralized algorithms, that can scale to very large networks. 
Dual decomposition was first applied to the NUM problem in 
[5], and has led to an extensive body of research on distributed 
algorithms for network optimization [6]-[8] and new ways to 
interpret existing network protocols [9]. 

Recent work by Zymnis et al. presented a specialized 
primal-dual interior-point method for the NUM problem [10]. 
Each Newton step is computed using the preconditioned 
conjugate gradient method (PCG). This proposed method had 
a significant performance improvement over the dual decom- 
position approach, especially when the network is congested. 



Furthermore, the method can handle utility functions which 
are not strictly concave. The main drawback of the primal-dual 
method is that it is centralized, while the dual decomposition 
methods are easily distributed. 

Other recent work by Bickson et al. [11] proposes an 
efficient way for computing the Newton step, which is the 
main computational effort of the primal-dual interior-point 
method using the Gaussian belief propagation (GaBP) algo- 
rithm, which is an efficient distributed algorithm. 

In the current paper we propose to combine both previous 
approaches. We present an efficient primal-dual interior point 
method, where the Newton step computed in each iteration is 
computed using the GaBP algorithm. Using extensive simula- 
tions with very large scale networks we compare the perfor- 
mance of our novel method to previous approaches including 
an interior-point method using PCG, and dual decomposition 
methods. Despite of being distributed, our new construction 
exhibits significant performance improvements over previous 
approaches. 

Furthermore, we provide the first comparison of perfor- 
mance of the GaBP algorithm vs. the PCG method. The 
PCG method is a state-of-the-art method used extensively in 
large-scale optimization applications. Examples include £i- 
regularized logistic regression [12], gate sizing [13], and slack 
allocation [14]. Empirically, the GaBP algorithm is immune to 
numerical problems with typically occur in the PCG method, 
while demonstrating a faster convergence. The only previous 
work comparing the performance of GaBP vs. PCG we are 
aware of is [15], which used a small example of 25 nodes, 
and the work of [16] which used a grid of 25 x 25 nodes. 

We believe that our approach is general and not limited to 
the NUM problem. It could potentially be used for the solution 
of other large scale distributed optimization problems. 

This paper is organized as follows. Section |ll] briefly 
overviews the NUM problem formulation. Section Hill outlines 
previous algorithms for solving the NUM problem, including 
dual descent and truncated Newton method. We present our 
new construction which utilizes the GaBP algorithm in Section 
HVl Section [V] provides simulation results comparing the 
performance of the GaBP based algorithm with the previous 
approaches. We conclude in Section [VTl 



II. Problem formulation 

There are n flows in a network, each of which is associated 
with a fixed route, i.e., some subset of m finks. Each flow 
has a nonnegative rate, which we denote /i, . . . , /„. With the 
flow j we associate a utility function Uj : R ^ R, which is 
concave and twice differentiable, with domUj C R^. The 
utility derived by a flow rate fj is given by Uj{fj). Tfie total 
utility associated with all the flows is then U{f) = C/i(/i) + 

■■■ + UMn). 

The total traffic on a link in the network is the sum of the 
rates of all flows that utilize that link. We can express the 
link traffic compactly using the routing or link-route matrix 
R e R"^", defined as 



Ri 



1 flow j's route passes over link i 
otherwise. 



Each link in the network has a (positive) capacity ci, . . . , c^. 
The traffic on a link cannot exceed its capacity, i.e., we have 
Rf < c, where < is used for componentwise inequality. 

The NUM problem is to choose the rates to maximize 
total utility, subject to the link capacity and the nonnegativity 
constraints: 



maximize U{f) 

subject to Rf < c, f > 0, 



(1) 



with variable / G R". This is a convex optimization problem 
and can be solved by a variety of methods. We say that / is 
primal feasible if it satisfies Rf < c, / > 0. 
The dual of problem ([T]) is 



minimize A^c + X;"=^(~i7j)*(-rjA) 
subject to A > 0, 



(2) 



where A G R™ is the dual variable associated with the capacity 
constraint of problem ([T]), rj is the jth column of R and 
{—UjY is the conjugate of the negative jth utility function 
[17, §3.3], 

{—Uj)*{a) = sup(ax + Uj{x)). 
We say that A is dual feasible if it satisfies A > and A G 

n;udom(-f/,)*. 

III. Previous work 

In this section we give a brief overview of the dual- 
decomposition method and the primal-dual interior point 
method proposed in [10]. 

A. Dual decomposition 

Dual decomposition [3]-[6] is a projected (sub)gradient 
algorithm for solving problem (|2|i, in the case when all utility 
functions are strictly concave. We start with any positive A, 
and repeatedly carry out the update 



A := 



argmax {Uj{x) - x{rj \)) , j = 1, . 

a;>0 



where a > is the step size, and x+ denotes the entrywise 
nonnegative part of the vector x. It can be shown that for small 
enough a, f and A will converge to and A*, respectively, 
provided all Uj are differentiable and strictly concave. The 
term s — c — Rf appearing in the update is the slack in the 
link capacity constraints (and can have negative entries during 
the algorithm execution). It can be shown that the slack is 
exactly the gradient of the dual objective function. 

Dual decomposition is a distributed algorithm. Each flow is 
updated based on information obtained from the links it passes 
over, and each link dual variable is updated based only on the 
flows that pass over it. 

B. Primal-dual interior point method 

The primal-dual interior-point method is based on using 
a Newton step, applied to a suitably modified form of the 
optimality conditions. The modification is parametrized by a 
parameter t, which is adjusted during the algorithm based 
on progress, as measured by the actual duality gap (if it is 
available) or a surrogate duality gap (when the actual duality 
gap is not available). 

We first describe the search direction. We modify the 
complementary slackness conditions to obtain the modified 
optimality conditions 

-VU{f) + R^X- fi = 

diag(A)s - (l/t)l 
diag(//)/ = (l/t)l, 

where t > is a parameter that sets the accuracy of 
the approximation. (As i ^ cx), we recover the optimality 
conditions for the NUM problem.) Here we implicitly assume 
that /, s, A, /i > 0. The modified optimality conditions can be 
compactly written as rt{f, A, /i) = 0, where 



rt{f,X,lJ.) = 



-VU{f) + R:^X- fi 
diag(A)s - (1A)1 
diag(Ai)/ - (1A)1 



The primal-dual search direction is the Newton step for 
solving the nonlinear equations rt(/, A,/Lt) = 0. If y — 
(/, A, fi) denotes the current point, the Newton step Ay — 
(A/, AA, A/i) is characterized by the linear equations 



niy + Ay) w rt(y) + r[{y)Ay 
which, written out in more detail, are 



0, 



-V^(7(/) 
^diag(A)i? 
diag(M) 



diag(s) 




-/ 


diag(/) 





" A/ " 




AA 




Afj. 



(A-a(c-i?/)), 



~rt{f,X,fJ.). 
(3) 

During the algorithm, the parameter t is increased, as the 
primal and dual variables approach optimality. When we have 
easy access to a dual feasible point during the algorithm, we 
can make use of the exact duality gap rj to set the value of t; 
in other cases, we can use the surrogate duality gap i). 

The primal-dual interior point algorithm is given in [17, 
§11.7], [18]. 



The most expensive part of computing the primal-dual 
search direction is solving equation Q. For problems of 
modest size, i.e., with m and n no more than 10'*, it can 
be solved using direct methods such as a sparse Cholesky 
decomposition. 

For larger problem instances [10] proposes to solve ^ ap- 
proximately, using a preconditioned conjugate gradient (PCG) 
algorithm [19, §6.6], [20, chap. 2], [21, chap. 5]. When an 
iterative method is used to approximately solve a Newton 
system, the algorithm is referred to as an inexact, iterative, 
or approximate Newton method (see [20, chap. 6] and its 
references). When an iterative method is used inside a primal- 
dual interior-point method, the overall algorithm is called 
a truncated-Newton primal-dual interior-point method. For 
details of the PCG algorithm, we refer the reader to the 
references cited above. Each iteration requires multiplication 
of the matrix by a vector, and a few vector inner products. 

IV. Our new construction 

Previous work of Zymnis et aJ. [10] shows that when ap- 
plying the interior-point Newton method to the NUM problem, 
each Newton step involves a solution of Eq. |3] where the 
solution (A/, AA, A/i)^ is the Newton search direction. 

Recent results by Bickson etal. [22], [23] utilizes the GaBP 
algorithm as an efficient distributed algorithm for solving a 
system of linear equations. For utilizing the GaBP algorithm, 
we first normalize Eq. ^ by (1. — 1/A, — l//i) to get the 
following equivalent system of linear equations: 



The corresponding graph of the covariance matrix A is Q, 
with edge potentials ('compatibility functions') ipij and self- 
potentials ('evidence') 1/)^. These graph potentials are deter- 
mined according to the following pairwise factorization of the 
Gaussian distribution p{x) oc nr=i 4>iixi) Hfij} '^iji^i^Xj), 
resulting in ^jjij{xi,Xj) = exp{~'XiAijXj), and (j)i{xi) = 
exp (fejXi — Aiixf/i). The set of edges corresponds 
to the set of non-zero entries in A (Eq. |5]l. Hence, we would 
like to calculate the marginal densities, which must also be 
Gaussian, 

p{xi)^M{^x^ = {A-H},,Pr^ = {A-^}u), 

where and Pi are the marginal mean and inverse variance 
(a.k.a. precision), respectively. Recall that, according to [24], 
the inferred mean /i is identical to the desired solution of (Eq. 
nil. The GaBP update rules are summarized in Table U We 
use the notation as the set of node i graph neighbors, 
excluding i. 







-J 




" A/ " 


R 


- diag(s/A) 







AA 


-I 





diag(//At) 




A^i 



^-rt{f,X, 
(4) 

where ri(/,A,/i) = ((-VC/(/) + - /i)^, (-s + 
(A/t))^, (-/+ (^/i))^)^. Note that the new system of linear 
equations is symmetric, a condition required by the GaBP 
algorithm. 

The formulation ^ allows us to shift the linear system of 
equations from an algebraic to a probabilistic domain. Instead 
of solving a deterministic vector-matrix linear equation, we 
now solve an inference problem in a graphical model describ- 
ing a certain Gaussian distribution function. Following [24] 
we define the joint covariance matrix 



# 


Stage 


Operation 


1. 


Initialize 


Compute Pa = Au and fiu = bi/Aa. 
Set Pki = and fiki = 0, Vfc £ 


2. 


Iterate 


Propagate Pki and pki, Vfc £ N(i) . 
Compute Pi\j = Pa + Y.ketm\j 

Compute Pij = —AijPT^.Aji 

Hij = ~Pij Aij^i\j. 


3. 


Check 


If Pij and Pij did not converge, 
return to #2. Else, continue to #4. 


4. 


Infer 


Pi ~ Pa + '}2keN(i) -Pfe* 

Pi ~ Pi {Piipii + X]fcgN(i) Pkipki) ■ 


5. 


Output 


Xi — Pi 



A^ 



R 
-I 



R^ 
-diag(s/A) 






diag(//M) 



(5) 



and the shift vector b ^ {{~VU{f) + - n)^ , {-s + 
{\/t)Y, (-/ + {n/t)YY. We further denote the search 
direction x = (A/^, AA^, A/x^)^. 

Given the covariance matrix A and the shift vector b, one 
can write explicitly the Gaussian density function 

p{x) ~ exp(-l/2x'^Ax + b^x) 

Now, we are interested in computing the MAP assignment: 

x* = argmaxexp(— l/2x^Ax + b^x) 



TABLE I: Computing x = A'^b via GaBP [23]. 

It is known that if GaBP converges, it results in exact 
inference [25]. Determining the exact region of convergence 
remain open research problems. All that is known is a suffi- 
cient (but not necessary) condition stating that GaBP converges 
when the spectral radius satisfies p{\Ik — A\) < 1 [26], [27]. 
A stricter sufficient condition [25], determines that the matrix 
A must be diagonally dominant (i.e., \Aii\ > X^j^^i l^jjli^*) 
in order for GaBP to converge. Recently, a new technique 
for forcing convergence for any column-dependent matrices is 
proposed in [28]. An upper bound on convergence speed is 
given in [11]. 

V. Experimental results 

A. Small experiment 

In our first example we look at the performance of our 
method on a small network. The utility functions are all 
logarithmic, i.e., Uj{fj) — log/j. There are n = 10'^ flows, 
and TO = 2 10'^ links. The elements of R are chosen randomly 
and independently, so that the average route length is 10 links. 
The link capacities Ci are chosen independently from a uniform 



distribution on [0.1, 1]. For this particular example, there are 
about 10^ nonzero elements in R (0.5% density). 

We compare three different algorithms for solving the NUM 
problem: The dual-decomposition method, a truncated Newton 
method via PCG and a customized Newton method via the 
GaBP solver. Out of the examined algorithms, the Newton 
method is centralized, while the dual-decomposition and GaBP 
solver are distributed algorithms. The source code of our 
Matlab simulation is available on [29]. 




Customized Newton via GaBP 
Trancated Newton via PCG 
Dual decomposition 



50 100 150 200 250 300 

Iteration number 



System size m=20,000 n=1 0,000 




Number of iterations 



Fig. 2: Convergence rate in the larger settings. 

10^** where the truncated Newton method implemented via 
PCG converged in 13 steps to the same accuracy. However, 
when examining the iteration count in each Newton step (the 
Y-axis) we see that the GaBP remained constant, while the 
PCG iterations significantly increase as we are getting closer 
to the optimal point. 



Fig. 1: Convergence rate using the small settings. 

Figure [T] depicts the solution quality, where the X-axis 
represents the number of algorithm iterations, and the Y-axis is 
the surrogate duality gap (using a logarithmic scale). As clearly 
shown, the GaBP algorithm has a comparable performance to 
the sparse Cholesky decomposition, while it is a distributed 
algorithm. The dual decomposition method has much slower 
convergence. 

B. Larger experiment 

Our second example is too large to be solved using the 
primal-dual interior-point method with direct search direction 
computation, but is readily handled by the truncated-Newton 
primal-dual algorithm using PCG, the dual decomposition 
method and the customized Newton method via GaBP. The 
utility functions are all logarithmic: Uj{fj) = log/j. There are 
n = 10* flows, and m — 2-10^ links. The elements of R and c 
are chosen as for the small example. For dual decomposition, 
we initialized all as 1. For the interior-point method, we 
initialized all A; and /i^ as 1. We initialize all fj as 7, where 
we choose 7 so that Rf < 0.9c. 

Our experimental results shows, that as the system size 
grows larger, the GaBP solver has favorable performance. 
Figure |2] plots the duality gap of both algorithms, vs. the 
number of iterations performed. 

Figure |3] shows that in terms of Newton steps, both methods 
had comparable performance. The Newton method via the 
GaBP algorithm converged in 11 steps, to an accuracy of 
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Fig. 3: Iteration count per Newton step. 

We have experimented with larger settings, up to 71 = 10^ 
flows, and m = 2 ■ 10^ links. The GaBP algorithm converged 
in 11 Newton steps with 7-9 inner iteration in each Newton 
step. The PCG method converged in 16 Newton steps with an 
average of 45 inner iterations. 

C. Numerical issues 

Overall, we have observed three types of numerical prob- 
lems with the PCG method. First, the PCG Matlab implemen- 
tation runs into numerical problems and failed to compute 



the search direction. Second, the line search failed, which 
means that no progress is possible in the computed direction 
without violating the problem constraints. Third, when getting 
close to the optimal solution, the number of PCG iterations 
significantly increases. 

The numerical problems of the PCG algorithm are well 
known, see of example [30], |31|. In contrary, the GaBP 
algorithm did not suffer from the above numerical problems. 

Furthermore, the PCG is harder to distribute, since in each 
PCG iteration a vector dot product and a matrix product 
are performed. Those operations are global, unlike the GaBP 
which exploits the sparseness of the input matrix. 

VI. Conclusion 

We propose an efficient distributed solution of the NUM 
problem using a customized Newton method, implemented 
via the GaBP algorithm. We compare the customized Newton 
method performance with state-of-the-art algorithms, includ- 
ing a dual descent method and a truncated Newton method, 
over large scale settings. We observe both faster convergence 
of the GaBP algorithm compared to both the preconditioned 
conjugate gradient and sparse Cholesky factorization. Further- 
more, the GaBP does not suffer from numerical problems 
which affect the performance of the preconditioned conjugate 
gradient method. 

We believe that the NUM problem serves as a case study 
for demonstrating the superior performance of the GaBP 
algorithm in solving sparse systems of linear equations. Since 
the problem of solving a system of linear equations is a 
fundamental problem in computer science and engineering, we 
envision many other appUcations for our proposed method. 
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